function S = myprepsimulate1(This,S,Opt)
% myprepsimulate1  [Not a public function] Prepare loop-independent data.
%
% Backed IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2015 IRIS Solutions Team.

%--------------------------------------------------------------------------

ny = length(This.solutionid{1});
nxx = length(This.solutionid{2});

% Bkw compatibility.
if isfield(Opt,'nonlinear') && ~isempty(Opt.nonlinear)
    % ##### Feb 2015 OBSOLETE and scheduled for removal.
    utils.warning('obsolete', ...
        ['Option nonlinear= is obsolete, and ', ...
        'will be removed from IRIS in a future release. ', ...
        'Use new options method= and nonlinper= instead.']);
    Opt.method = 'selective';
    Opt.nonlinper = Opt.nonlinear;
end

S.Method = lower(Opt.method);
S.Solver = lower(Opt.solver);
S.Display = lower(Opt.display);
S.IsDeviation = Opt.deviation;
S.IsAddSstate = Opt.addsstate;
S.IsError = Opt.error;
S.Eqtn = This.eqtn;
S.NameType = This.nametype;
S.IxYLog = This.IxLog(real(This.solutionid{1})).';
S.IxXLog = This.IxLog(real(This.solutionid{2})).';

S.NPerNonlin = 0;
if isequal(S.Method,'selective')
    if any(This.IxNonlin)
        S.NPerNonlin = Opt.nonlinper;
        if isequal(S.NPerNonlin,@all)
            S.NPerNonlin = S.NPer;
        end
        S.ZerothSegment = 0;
        if S.NPerNonlin==0
            S.Method = 'firstorder';
        end
    else
        utils.warning('model:myprepsimulate1', ...
            ['No nonlinear equations marked in the model file. ', ...
            'Switching to first-order simulation.']);
        S.Method = 'firstorder';
    end
end

if isequal(S.Method,'selective')
    S.Selective = struct();
    S.Selective.IxNonlin = This.IxNonlin;
    S.Selective.Tolerance = Opt.tolerance;
    S.Selective.MaxIter = Opt.maxiter;
    S.Selective.Lambda = Opt.lambda;
    S.Selective.UpperBnd = Opt.upperbound;
    S.Selective.IsFillOut = Opt.fillout;
    S.Selective.ReduceLmb = Opt.reducelambda;
    S.Selective.MaxNumelJv = Opt.maxnumeljv;
    label = This.eqtnlabel;
    ixEmpty = cellfun(@isempty,label);
    label(ixEmpty) = This.eqtn(ixEmpty);
    S.Selective.EqtnLabelN = label(This.IxNonlin);
    [S.Selective.EqtnN,S.Selective.EqtnNI] = doEqtnN(This.IxNonlin);
    if isequal(S.Solver,@auto) ...
            || ( ischar(S.Solver) && strcmpi(S.Solver,'qad') )
        S.Solver = @qad;
    end
    S.Display = S.Display;
    if isequal(S.Solver,@qad)
        if isequal(S.Display,'off') || isequal(S.Display,'on')        
            S.Display = 0;
        elseif isequal(S.Display,@auto) || ~isnumericscalar(S.Display)
            S.Display = 100;
        end
    else
        if isequal(S.Display,@auto)
            S.Display = 'iter';
        end
    end
    S.Selective.Display = S.Display;
    % Positions of variables in [y;xx] vector that occur in nonlinear
    % equations. These will be  combined with positions of exogenized variables
    % in each segment.
    ixXUpdN = false(nxx,1);
    for iEq = find(This.IxNonlin)
        tok = regexp(This.EqtnN{iEq},'\<xx\>\((\d+),','tokens');
        tok = [tok{:}];
        ix = eval(['[',sprintf('%s,',tok{:}),']']);
        ixXUpdN(ix) = true;
    end
    S.Selective.IxUpdN = [false(ny,1);ixXUpdN];
    S.Selective.NOptimLambda = double(Opt.noptimlambda);
    S.Selective.NShanks = Opt.nshanks;
    
elseif isequal(S.Method,'global')
    S.NPerNonlin = S.NPer; % NPerNonlin is used in MYPREPSIMULATE2.
    S.Global = struct();
    ix = This.eqtntype==1 | This.eqtntype==2;
    [S.Global.EqtnN,S.Global.EqtnNI] = doEqtnN(ix);
    S.Global.IxXCurr = ( imag(This.solutionid{2})==0 ).';
    if isequal(S.Solver,@auto)
        S.Solver = @fsolve;
    end
    if isequal(S.Display,@auto)
        S.Display = 'final';
    end
    S.Global.Unlog = Opt.unlog;
    S.Global.Tolerance = model.DEFAULT_STEADY_TOLERANCE;
    S.Global.Times = Opt.times;
    if isequal(S.Global.Times,@auto)
        S.Global.Times = 10*(ny+nxx);
    end
end

% Steady state update.
S.IsUpdate = Opt.sstateupdate ...
    && ~isempty(This.PosUpdate) && any(This.PosUpdate>0);
if S.IsUpdate
    % Initialize and preprocess sstate, chksstate, solve options.
    S.Update = struct();
    S.Update.sstate = mysstateopt(This,'silent',Opt.sstate);
    S.Update.chksstate = mychksstateopt(This,'silent',Opt.chksstate);
    S.Update.solve = mysolveopt(This,'silent',Opt.solve);
    ixUpdate = This.PosUpdate>0;
    eqtnUpd = This.EqtnN(This.eqtntype==5);
    eqtnUpd = eqtnUpd(ixUpdate);
    eqtnUpd = [ model.EQTN_NONLIN_PREAMBLE, '[', eqtnUpd{:}, ']' ];
    eqtnUpd = str2func(eqtnUpd);
    S.Update.EqtnN = eqtnUpd;
    S.Update.PosUpdate = This.PosUpdate(ixUpdate);
end

if isequal(S.Solver,@lsqnonlin) || isequal(S.Solver,@fsolve)
    Opt.solver = S.Solver;
    Opt.display = S.Display;
    Opt.tolfun = model.DEFAULT_STEADY_TOLERANCE;
    Opt.tolx = model.DEFAULT_STEADY_TOLERANCE;
    [~,S.OptimSet] = irisoptim.myoptimopts(Opt);
end

return



    function [EqtnN,EqtnNI] = doEqtnN(Ix)
        EqtnN = This.EqtnN(Ix);
        EqtnN = [ model.EQTN_NONLIN_PREAMBLE, '[', EqtnN{:}, ']' ];
        if true % ##### MOSW
            EqtnN = str2func(EqtnN);
        else
            eqtnN = mosw.str2func(eqtnN); %#ok<UNRCH>
        end
        % Individual equations to catch errors.
        EqtnNI = cell(size(This.eqtn));
        for ii = find(Ix)
            EqtnNI{ii} = This.EqtnN{ii};
            EqtnNI{ii} = [ model.EQTN_NONLIN_PREAMBLE, EqtnNI{ii} ];
            if true % ##### MOSW
                EqtnNI{ii} = str2func(EqtnNI{ii});
            else
                eqtnNI{ii} = mosw.str2func(eqtnNI{ii}); %#ok<UNRCH>
            end
        end
    end % doEqtnN()
end
